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The Wigner d function, which is the essential part of an irreducible representation of SU(2) and 
SO(3) parameterized with Euler angles, has been know to suffer from a serious numerical errors at 
high spins, if it is calculated by means of the Wigner formula as a polynomial of cos and sin of 
half of the second Euler angle. This paper shows a way to avoid this problem by expressing the d 
functions as the Fourier series of the half angle. A precise numerical table of the coefficients of the 
series is provided as Supplemental Material. 


I. INTRODUCTION 


where nmin and n„iax are zero or positive integers. 


The matrix elements of the rotation operator between 
angular-momentum eigenstates are called the Wigner D 
function [I| . When the rotation is specified with the three 
Euler angles (</>, 9, ip), and the eigenstates are labeled 
with the magnitude (j = 0, 1, • •) and the z com¬ 

ponent (m or k = —j, — J + 1, • • •, j) of the angular mo¬ 
mentum vector, the D function can be decomposed into 
three factors, 

= ( 1 ) 

where 

( 2 ) 

is the nontrivial part which needs some method to eval¬ 
uate. It is called the Wigner (small) d function. In the 
standard phase convention for angular-momentum eigen¬ 
states, the matrix elements of jy are purely imaginary 
and thus takes on real numbers. 

One may be anxious about the fact that Bohr and 
Mottelson[3 define the rotation matrix as the 

complex conjugate of the right-hand side (r.h.s.) of 
Eq.([T]). However, their definition for is identical with 
Wigner’s and the definition of the d function is unique. 

The Wigner D function is used in various fields of 
physics. In some applications, those for large values of j 
(say, > 50) are necessary. An example in nuclear struc¬ 
ture physics is the projection from spatially deformed 
solutions of modern realistic mean-field models to eigen¬ 
states of large angular momentum. For toy models, too, 
one occasionally needs d functions for very large j to 
confirm the validity of one’s picture under extreme con¬ 
ditions. (e.g., a two-rotor model of Refs. Q and i)- 
The explicit form of the d function is given by the 
Wigner formula[lj, which can be written as, 

'^max 
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nmin = max(0, k — m), (4) 

«max = min(j - m,j-I-fc), (5) 


and 


Wr%6)=w: 
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( 6 ) 


jmk ^ VU + ^)Kj - m)Kj + - k)l , . 

{j — m — n)\{j -\- k — n)l{n -\- m — k)lnl 

However, this formula suffers from a serious loss of pre¬ 
cision at high spins (i.e., for large j) except in the neigh¬ 
borhood of 9 = 0,7r. 

For example, assuming that j is a positive integer, 9 = 
and m = k = 0, one obtains 



( 8 ) 


which, if j is even, becomes maximum at n = ^, 



1 

2^' [ij/2y.f 



(9) 


(The Stirling’s formula is used in the last approximation.) 
The absolute value of the d function is not greater than 
one because it is a matrix element of a unitary operator 
between normalized states. Wigner’s formula expresses 
the d function as a result of cancellation among terms 
of possibly huge size, ^ 2P For j ^ 54, the pre¬ 

cision of double-precision floating-point numbers (53-bit 
mantissa) is lost completely, and even quadruple preci¬ 
sion float numbers (113-bit mantissa) is lost completely 
for j ^ 114. 

A few remedies have been proposed [1, Q but the re¬ 
sults are not completely free of the precision loss. In this 
paper, I investigate the details of this loss of significance, 
and then present a perfect remedy to avoid such numer¬ 
ical difficulty. 
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II. THE FOURIER-SERIES EXPRESSION OF d 
FUNCTIONS 


half-integer values), one can express the coefficients 
by an integral 


One can see easily that terms like cos''' | sin^ | appear¬ 
ing in the r.h.s. of Eq. where A and p, are zero or pos¬ 
itive integers such that X + fJ, < 2 j, can be expressed as 
linear combinations of terms sin ^ (when fj, is odd) and 
cos ^ (when /i is even) with integers k in 0<K<X + fj, 
and K = X + (mod 2), by means of repeated applica¬ 
tions of the elementary trigonometric identities called the 
product-to-sum identities or the prosthaphaeresis formu¬ 
las. 

Because the power of sin ^ is fj, = m— k + 2n in Eq. ©, 
one can see that fj, = m— k (mod 2) and that the d func¬ 
tion is an even (odd) function if m—k is even (odd), which 
can be expanded only with cos (sin) function. This may 
also be deduced from Eq. (IA3I) . one of the properties of 
the d function which I have enumerated in the appendix 

E 

From these considerations, one can conclude that the 
Fourier expansion of the d function has the following 
form: 


fjmk 


1 

27r(l -I- 6,yo) 



dl,{9)f{u0)d0, 


(15) 


where 6,^0 = 1 for v — 0 and <5,^0 = 0 for = A, 1, |, • • •. 

By substituting the d function in Eq. m with Eqs. (H 
ED, I have derived a more useful expression for con¬ 
taining only four elementary operations of arithmetic. 
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(16) 


where nmin and Umax are those already defined by 
Eqs. @ and ®, the square brackets are the floor func¬ 
tion, i.e., [I + x] = I for integer I and real x in [0,1), 


p=\m — k\ (mod 2), (17) 


i.e.. 


= (10) 

where 

/={s}'»M 

and the summation runs over 


p = 0 for fc = TO, m ± 2, TO ± 4, • • •, (18) 

p=l for fc = TO ± 1, TO ± 3, • • •, (19) 

and 

p27Z 

Ixfj. = / cos"''a: sin^ X dx (20) 

Jo 

with zero or positive integers for A and p. If both A and 
p are even. 


k* — k'min: ^min “t“ 1? * * * ; j) (1^) 

with the values of t'min given in Table ID 

TABLE 1. The values of I'min appearing in Eq. GQl) 



Even m — k 

Odd m — k 

Even 2j 

0 

1 

Odd 2j 

1 

_2_ 

1 

_2_ 


For example, for J = | and to = — fc = 1, the Wigner 
formula ([3]) gives an expression, 

dl _i{9) = sin^ I — 12 cos^ | sin^ f + 18 cos"^ | sin^ | 


—4 cos® I sin |, (13) 

which can be rewritten in the form as 

1 , , 35 sin ^ — 5 sin ^ -I- 15 sin ^ — 9 sin | 

dl _.(0) =-^2 

2 ' 2 64 


(14) 

By utilizing the orthogonality of cos iy0 and sin vd over 
0 < d < 47r (considering that z/ can take both integer and 


-kxfi — 


27r(A- l)!!(Ai- 1)!! 

^+7on 


( 21 ) 


while I\i^ = 0 otherwise. 

Unlike the r.h.s. of Eq. ®, where the terms can have 
huge sizes and thus the numerical error is a serious prob¬ 
lem, the r.h.s. of Eq. (nni) is a summation of terms of 
order one or less and hence the problem is expected to 
disappear. This can be seen by calculating the integrals 
of the squares of the both sides of Eq. (El, 

/ di^,{0fd0 = Y,Y.*i"''*r' f{^0)f{d0)d0 

Jo u t, Jo 

= (22) 

Because the absolute values of d functions are < 1, the 
left-hand side is < 47r and, consequently, it holds < 

2 — S^Q. 

A further study from the numerical point of view has 
indicated that the maximum (among all the possible com¬ 
binations of TO, fc, and ix) value of is 1 for j < 1, 

decreases as j increases from an integer to the next half 
integer, and does not change as j increases from a half 
integer to the next integer. For the interval 50 < j < 100, 
the maximum value for integer j behaves as ~ 1.13/y^. 
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III. COMPUTATION OF THE NUMERICAL 
VALUES OF THE COEFFICIENTS 


Unfortunately, Eq. (HH) also suffers from a serious loss 
of significant digits in ordinary floating-point numerical 
calculations. Indeed, for even integer j, the term having 
the maximum magnitude among those in the r.h.s. of 
Eq. (ITSl) occurs at m = fc = 0, n = r = ij, v = j, with 
the maximum value of 


joo 

'^j/2 




220 + 1 ) 


(23) 


which is roughly 2^ times as large as the value given by 

Eq. (O. 

To avoid this problem, I first evaluated the r.h.s. of 
Eq. (1161) rigorously as rational numbers or square root 
of rational numbers by means of a formula-manipulation 
software MAXIMA. Numerical values to be used in pro¬ 
grams coded in Eortran, C, etc., can be calculated from 
such rigorous numbers to the full precision of the 64-bit 
floating-point number. 

However, the computation time turned out to be exces¬ 
sively long for large values of j. To speed up the computa¬ 
tion, I have changed the method to evaluate Eq. (HU) not 
rigorously but in terms of high-precision floating-point 
numbers (of MAXIMA). This does not seem to be a ma¬ 
jor drawback because numerical values are sufficient for 
most of practical purposes. 

I change the precision of floating-point numbers de¬ 
pending on j in such a way that the number of digits 
equals the common logarithm of the value of Eq. (ED 
divided by 10“^®. It increases with j, reaching 74 digits 
for j = 100. Precise 64-bit floating-point numbers can be 
obtained simply by truncating the high-precision results. 

Empirically, the time necessary to compute all the co¬ 
efficients for each j increases as j'^. The new method 
takes 43 h for j = 100 with a personal computer with a 
CPU Intel core-i7 3960X running at 3.3GHz, using one 
physical core. 

I use the obtained 64-bit floating-point number coeffi¬ 
cients to evaluate the Fourier-series formula for d func¬ 
tions. I provide data files of the numerical values of the 
coefficients, together with a sample FORTRAN90 pro¬ 
gram to read the data and calculate the values of the d 
function, as Supplemental Material Q- 

For each value of j, there are {2j -|-1)^ possible combi¬ 
nations of the values of m and k (because —j < m < j, 
~j ^ k < j). Only about a quarter of them are in¬ 
dependent, however, because of the properties of the d 
function expressed by Eqs. (IA4| . (IA6|) . and dm in the 
appendix E Therefore, I consider only such combina¬ 
tions as m > 0 and k < \m\ in the following analysis of 
the numerical precision. 

In other words, the coefficients have the following 
symmetry, 

^ ( 24 ) 

^ ( 25 ) 


^ (26) 

Hence, the numerical data for in the Supplemental 
Material are given only for m > 0 and k < \m\. 

There are ~ ^(jmax + |)^ coefficients for j = 
0, 5 ,1, • • •, jmax- The size of the memory to store them as 
64-bit floating-point numbers amounts to 27 (404) MiB 
for jmax=50 (100). My data are given as text files for the 
sake of compatibility, whose sizes are not very different 
from the above memory sizes after they are compressed 
(with the software GZIP). 

The magnitude of the coefficient becomes smaller 
for larger values of j. However, even at j = 100, 80% 
(90%) of the coefficients are larger than 10“® (10“^°). 

Some of the coefficients vanish according to rules un¬ 
mentioned so far. For example, = 0 if j is an integer, 
TO = 0 and/or k = 0, and j — v = 1 (mod 2). One can 
prove this vanishment by using Eq. (IA8I) . I do not use 
these additional rules but simply give zero values in the 
data files. 

The evaluations of cos and sin i/9 should be calcu¬ 
lated by means of a recursion relation, 

( cos {v + 1)9 ]_ I cos 9, —sin 9 \ I cos 1^9 \ 
sin(j^ -I- 1)9 J y sin 9, cos 9 y y sin y 

(27) 

which is nothing but the trigonometric (angle) addition 
theorem. For integer j, the initial values are cos 0 = 1 
and sin 0 = 0. For half integer j, one has to calculate, 
first, the initial values cos | and sin | and, second, cos 9 
and sin 9 using identities cos 9 = cos^ | — sin^ | and 
sin d = 2 sin | cos |. 



FIG. 1. Maximum and root-mean-square errors in the val¬ 
ues of cos and sinuO, numerically calculated with 64-bit 
floating-point numbers, over 0 < 8 < tv. The values of 9 are 
sampled with a uniform spacing of 10“® degree. The abscissa 
is = 0, |, 1, • • •, 100. Recursion means the usage of Eq. 
together with the initial values described in the following sen¬ 
tences. 

Straightforward evaluation of cos v9 and sin v9 (i.e., 
passing the value of v9 to the internal functions cos and 
sin) requires about j calls to the functions and is compu¬ 
tationally very inefficient. Moreover, as shown in Fig. [TJ 
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such straightforward evaluation causes slightly larger nu¬ 
merical errors probably due to the loss of significant digits 
in reducing the value of v9 to an interval such as [—f, f ], 
especially when the value of \ v9\ is large. The reason why 
the recursion formula (l27l) does not suffer from large er¬ 
rors even after hundred steps may be attributed to the 
fact that the magnitudes of cos and sin functions are al¬ 
ways not greater than one. 


IV. PRECISION OF THE NUMERICAL 
VALUES OF THE d FUNCTION 
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FIG. 2. The common logarithm of the error in the 64-bit 
floating-point numerical value of the d function d^^f,{9) as a 
function of m and k for j = 40 and 0 = 30°. The value of the 
d function is calculated by means of the Wigner formula in 
part (a) and the Fourier series expression in part (b). Only a 
quarter of the possible combinations of m and k are plotted 
because the error is symmetric in the lines m = ±fc. Errors 
smaller than 10“^° are painted with the same color as that 
for the error of 10“^°. 
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FIG. 3. The same as in Fig. [2]but for 9 = 60°. 
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FIG. 4. The same as in Fig. [2]but for 9 = 90°. 


In this section, I compare the errors of the values of 
calculated with 64-bit floating-point numbers ac¬ 
cording to the Wigner formula ([3]) and the Fourier series 
expression m- The errors have been calculated as the 
differences from the exact values calculated by applying a 
formula-manipulation software MAXIMA to the Wigner 
formula. 

Figures m [31 and m show the errors for 9 = 30°, 60°, 
and 90°, respectively. The errors are expressed as a func¬ 
tion of {m,k) while j is fixed at 40. I have found that 
similar plots for 9 > 90° look almost indistinguishable 
from those for 180° — 9 except that the sign of k is re¬ 
versed, as could be foreseen from Eq. (1X91) . 

One can see that the Wigner formula results in very 
large errors, up to ~ 10“® for 9 ~ 90° and m ^ k ^ 0, 
while the Fourier series expression gives precisely 15 dig¬ 
its irrespective of the values of 9, j, m, and k. I confi¬ 


dently recommend the Fourier series expression over the 
Winger formula already at j ^ 40. 

For a special purpose, however, the Wigner formula 
still has an advantage. For regions of the arguments 9 ^ 
0° (180°) and \m + k\ ^ 0 (2j), the Wigner formula 
has smaller errors than 10“^®, i.e., than the level of the 
almost constant error of the Fourier series expression. In 
such regions of the arguments, the magnitude of d?^f.{9) is 
very small and a small number of terms dominate in the 
summation of the Wigner formula, while many terms of 
order 1 cancel among themselves to give the small value 
in the Fourier series expression. Therefore, if one needs 
such high precision at those regions, it would be a good 
idea to develop a program which switches between the 
two formulas depending on the values of 0, j, m, and, k. 

In Fig. El I compare the maximum errors of cP^j^i9) 
as functions of j. The maximum is taken over the val- 
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FIG. 5. The common logarithm of the maximum error in the 
numerical value of the d function versus j calculated 

with the Wigner formula (blue solid line) and the Fourier 
series expression (red dashed line). The maximum is taken 
over the values of m, k, and 9 for each j. 

ues of 9 from 0° to 180° with an increment of 5° and all 
possible combinations of m and k. For both formulas, 
the maximum error increases as exponential functions of 
j. The Wigner formula increases its error, however, far 
faster than the Fourier series expression. For the interval 
20 < j < 100, the error of the Wigner formula can be 
approximated by log^g |error| Ri 0.294j — 17.2 while that 
of the Fourier series expression can be approximated by 
logio |error| r; O.OOGj — 14.8. In other words, by increas¬ 
ing j by one, the error of the Wigner formula is doubled, 
while that of the Fourier series expression increases only 
by 1.4%. 


V. SUMMARY 

I have shown that the Wigner formula for the d func¬ 
tion results in intolerable large numerical errors for large 
values of the angular momentum quantum number j. On 
the other hand, the Fourier series expression for the d 
function is shown to be free of such errors, providing 
precision of ~ 10“^'* even at j = 100. An analytic ex¬ 
pression for the coefficients of the Fourier series is given. 


Their numerical values, which are precise as far as 64-bit 
floating-point numbers can express, are provided as elec¬ 
tric files in Supplemental Material. Sample programs in 
FORTRAN90 to use the data files are also provided. 

Appendix A: Properties of the d function utilized in 
this paper 

I enumerate the symmetries of the d function to be 
referred to in this paper. 

First, the unitarity of rotations (i.e, the Hermite con¬ 
jugate operator is the inverse operator) means 

dLki(^) = dU-(^)- _(A1) 

Second, the composition of two rotations can be rewritten 
as multiplication of matrices to represent them, 

dl.iei +92)= ^2 dLAdi)di,{92). (A2) 

'^=-3 


I need three more relations, whose easiest derivation may 
be to use the Wigner formula (l3|) as in Ref. [l|. 


dLkid) = 

(A3) 

dLkid) = dl^^_^{9), 

(A4) 

dLkM = 

(AS) 

From Eqs. (lAlIl and (jA.3ll. one can prove 

dUd) = 

(A6) 

from Eqs. (IA4I1 and (|A6p. 

dLkid) = 

(A7) 

and from Eqs. (IA2I1. (IA5I) . and (IA6I). 


(AS) 

dU^-d) = i-l)^^^d2-kid)- 

(A9) 
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